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Abstract 

Background: The human ether-a-go-go related gene 1 (hERGl), which codes for a potassium ion channel, is a key 
element in the cardiac delayed rectified potassium current, iKr, and plays an important role in the normal repolarization 
of the heart's action potential. Many approved drugs have been withdrawn from the market due to their prolongation 
of the QT interval. Most of these drugs have high potencies for their principal targets and are often irreplaceable, thus 
"rehabilitation" studies for decreasing their high hERGl blocking affinities, while keeping them active at the binding sites 
of their targets, have been proposed to enable these drugs to re-enter the market. 

Methods: In this proof-of-principle study, we focus on cisapride, a gastroprokinetic agent withdrawn from the market 
due to its high hERGl blocking affinity. Here we tested an a priori strategy to predict a compound's cardiotoxicity using 
de novo drug design with molecular docking and Molecular Dynamics (MD) simulations to generate a strategy for the 
rehabilitation of cisapride. 

Results: We focused on two key receptors, a target interaction with the (adenosine) receptor and an off-target interaction 
with hERGl channels. An analysis of the fragment interactions of cisapride at human A2A adenosine receptors and hERGl 
central cavities helped us to identify the key chemical groups responsible for the drug activity and hERGl blockade. A set 
of cisapride derivatives with reduced cardiotoxicity was then proposed using an in-silico two-tier approach. This set was 
compared against a large dataset of commercially available cisapride analogs and derivatives. 

Conclusions: An interaction decomposition of cisapride and cisapride derivatives allowed for the identification of key 
active scaffolds and functional groups that may be responsible for the unwanted blockade of hERGl. 

Keywords: hERG K channel, A2A adenosine receptor, 5HT-4 receptor. Molecular docking. Molecular dynamics simulations. 
Drug design. Drug databases 



Background 

Several classes of potassium channels are involved in regu- 
lating the heart rate by setting the amplitude and duration 
of the action potential and the resting membrane poten- 
tial Abnormalities in the function of these ion channels 
due to inherited mutations or pharmacological blockage 
can prolong the duration of the action potential, leading 
to the development of severe arrhythmias (i.e., long QT 
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syndromes - LQTS). Genetic analysis has revealed that 
mutations in potassium channels, such as the human 
ether- a-go-go related gene (hERG) and KvLQTl, establish 
a molecular basis for LQTS [1-4]. A growing number of 
diseases have also been linked to genetic mutations in 
potassium channels. The channelopathies related to po- 
tassium channels include various cancer types, type 2 
Bartter s syndrome, type 1 episodic ataxia, and hyper- 
insulinemic hypoglycemia [4]. The best-known feature 
of the hERGl channel is its unique promiscuity in bind- 
ing to a wide range of organic molecules. A broad panel 
of organic compounds used in common cardiac and 
non-cardiac medications (e.g., antibiotics, antihista- 
mines and antibacterial agents) are thought to cause a 
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reduction in the repolarizing current Ixr by blocking the 
central cavity of hERG and similar channels, leading 
to ventricular arrhythmia [5]. Several Food and Drug 
Administration (FDA) approved drugs (i.e., terfena- 
dine, cisapride, astemizole and grepafloxin) have been 
withdrawn from the market, while others like thiorida- 
zine, haloperidol, sertindole, and pimozide have been 
restricted due to their effect on the function of the 
hERG channel. The discovery of drug-related arrhyth- 
mias has led to mandatory drug screening for hERGl 
blockage by both the FDA and the European Medicines 
Agency (EMEA). Because most of these drugs have 
high binding affinity profiles for their principal targets, 
"rehabilitation" studies aimed at reducing their side 
effects (i.e., decreasing their high hERGl blocking af- 
finities) while keeping them efficient at binding to their 
original targets have become increasingly common. 
These studies may allow these drugs to reenter the 
market. 

We chose cisapride to study as a model hERG blocker 
that can potentially be rehabilitated and returned back to 
the market. Researchers at Janssen Pharmaceutica discov- 
ered cisapride (a gastroprokinetic agent; trade names: pre- 
pulsid, propulsid) in 1980. Cisapride, a prokinetic agent, 
increases gastrointestinal motility and acts as a selective 
serotonin agonist for 5HT-4 receptors. It also relieves 
gastrointestinal symptoms (i.e., constipation and bloating) 
by indirectly stimulating the release of acetylcholine in 
muscarinic receptors. In the gastrointestinal tract, the acti- 
vation of these receptors stimulates smooth muscle con- 
traction. While cisapride's main target is thought to be the 
5-hydroxytryptamine receptor 4 (5HT-4), it can selectively 



target a number of other G -protein-coupled receptors 
(GPCRs). For example, cisapride has binding affinities for 
5HT-4 (-14 nM) and adrenergic receptors (16 nM at the 
a-1 adrenergic receptor). Cisapride was marketed in the 
USA from 1993 to 2000; the use of cisapride has been 
associated with 341 reports of cardiovascular problems, 
including LQTS and torsades-de-pointes, and 80 re- 
ports of death [6] . After seven years on the market, cisa- 
pride was withdrawn in the USA and was limited for 
use in many other countries due to its high hERGl 
blocking affinity (IC50 6.5 nM, Mohammad et al [7]; 
IC50 44.5 nM, Rampe et al [8,9]; IC50 15.0 nM Drolet 
et al, [10]). A goal of our proof-of-principle study is to 
establish the feasibility of an in silico cardiotoxicity as- 
sessment with a multi-target computational approach 
for rehabilitation that: 

(i) Assesses drugs for their hERG -blocking ability; 

(ii) Identify the active components responsible for the 
original target activity; 

(iii) Caps or modifies the moieties responsible for the 
hERG blockade. 

The workflow chart is shown in Figure 1. We combined 
step-by-step ligand modification focusing on key func- 
tional groups (Figure 2), two-target receptor docking, 
and molecular dynamics (MD) simulations to achieve this 
goal. First, we determined the key-molecular fragments of 
cisapride responsible for its high-affinity binding to the 
A2A receptor using all-atom MD simulations and Mole- 
cular Mechanics/Generalized Born Surface Area (MM/ 
GBSA) binding energy decompositions, which are similar 



Determination of key-molecular Tragments of 
Cisapride Responsible for high-affinity binding to A2A 
Using all-atom MD simulations and MM/GBSA binding 
energy decompositions 



Micro-modifications of Cisapride at the specific 
regions that responsible for the hERG blockade and a 
less significant role in the drug binding at A?a receptor 



Docking of Cisapride derivatives at the hERG central 
cavity, compounds thai have low-blocking affinities 
are selected for further tests 



Screening of Cisapride-! ike synthesized compounds 
from ZINC data base at hERG central cavity and Aja 
receptor 



To investigate the similarities of binding pockets of 
A2A and 5HT4 GPCRs. homology modeling of 5-HT4 
receptor model is performed and predicted binding 
energies of known drugs are compared at both targets 



List of rehabiittated Cisapride derivatives 



Figure 1 The block-scheme for the work-flow in the computational rehabilitation of cisapride derivatives. 
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Figure 2 The group definition of cisapride (i.e., R1 = -C6H4-F; 
R2 = -O-CH2-; R3 = -CH2-CH2-; R4 = -C5H3N-O-CH3; R5 = -NH-CO-; 
R6 = C6H2-CI-NH2; R7 = -O-CH3) for per-group decomposition 
and inter-monomer contacts from the MM/GBSA analysis. 



to the approaches used in previous studies [11-13]. Next, 
small modifications of the sites responsible for the 
hERG blockade that played a less significant role in the 
stabilization of cisapride in the binding pocket of the A2A 
receptor were made. Our previously developed atomistic 
models of the open-state hERGl channel were used to 
guide the modifications of cisapride. The most potent 
drug variants were tested for their binding to the cen- 
tral cavity of the hERGl channel. The compounds that 
showed low affinities were then selected for further 
analysis. Because no crystal structure of the 5HT-4 re- 
ceptor is currently available, we used a recently crystal- 
lized agonist-bound active GPCR from the same family 
as the target (the human adenosine A2A receptor, PDB 
ID: 3QAK) [14]. Since cisapride can selectively bind to 
a number of other GPCRs (including A2A receptor) 
with a similar range of binding affinity at a similar 
binding pocket, we chose to work with the available 
agonist-bound active target site of the A2A receptor to 
screen cisapride and its derivatives. Prior binding site 
[15] and comparative homology modeling studies for 
these receptors indicate a high similarity between the 
5HT-4 and A2A binding pockets. Finally, a list of pos- 
sible modifications to the original cisapride molecule 
was generated. MD simulations with MM/GB-SA com- 
putations, database drug screening and de novo design 
studies clearly showed that the shorter alkyl chain in 
cisapride analogues are key to retaining their binding 
to the A2A receptor while remediating the blockade 
of the hERGl channel. To compare in silico results to 
de novo developments we screened a large panel of 
already synthesized cisapride analogs. Small molecule 
databanks (i.e., ZINC [16]) were screened for synthesized 
cisapride derivatives and the relevant literature was re- 
viewed to identify their activity in both the GPCR and 
hERG targets. Dual target docking combined with all- 
atom MD simulations were used to establish the key 



interactions between cisapride and its derivatives re- 
sponsible for their high-affinity binding to the A2A and 
how they receptor triggering hERGl blockade. Using 
this combination of techniques it is possible to assess 
novel compounds a priori for their cardiotoxicity risks 
associated with hERGl blockade as well as to identify 
sets of functional groups responsible for on-target 
binding. 

Methods 

Molecular docking 

The structures were initially optimized using Schrodinger s 
Macromodel module to perform an all-atom MM geom- 
etry optimization. The optimization was carried out with 
the OPLS 2005 force field and the Polak-Rebiere Conju- 
gate Gradient (PRCG) energy minimization method with a 
0.001 kcal mol"^ A"^ energy gradient convergence criter- 
ion [17]. The resultant ligand structures were docked to 
the targets using the following docking algorithms: Glide/ 
Induced Fit Docking (IFD) [18], FlexX [19], Autodock 
[20], and Generalized Optimized Ligand Docking (GOLD) 
[21]. The details of the docking algorithms are described 
below: 

Glide/IFD: The Glide-XP (extra precision) (v.5.0) and 
Induced Fit Docking (IFD) modules of the Maestro suite 
were used for the docking calculations. The docking 
studies were performed with the following steps: (i) 
constrained minimization of the receptor with an root 
mean square deviation (RMSD) with a cutoff of 0.18 A; 
(ii) initial Glide docking of each Ugand using soft po- 
tentials; (iii) refinement of the derived docking poses 
(i.e., minimization of the docking poses within 20 A of 
the ligand poses) with Schrodinger s Prime module; and 
(iv) Glide re-docking of the protein-ligand complexes. 
GOLD: The GOLD program (v.5.0. 1) was used with two 
default docking scores (GOLD Fitness and ChemScore). 
Part of the receptor's flexibility was accounted for by 
assigning flexibility to 10 specific amino acid residues at 
the active site. The side chains of these amino acid resi- 
dues were selected as flexible rotamers. The rotamers 
progressed by 10° increments to cover a full 360° rota- 
tion. The default genetic algorithm parameters (100 for 
the population size, 5 for the number of islands, 100000 
for the number of genetic operations and 2 for the niche 
size) were used. However, the maximum number of 
runs was set to 100 for each docking simulation. FlexX: 
The FlexX program (v.4.0) from BioSolvelT was also 
used. The default algorithm parameters were used for 
the docking and construction of the active sites of the 
receptor. The solutions per fragment and per iteration 
were both set to 2000. AutoDock (v.4.0): The number of 
grid points in each direction was 126 with a grid spacing 
of 0.4 A. The number of hybrid Genetic Algorithm- 
Local Search (GA-LS) runs was 200. 
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Homology modeling 

The SWISS-MODEL homology modeling program [22] 
was used for the development of the 5HT-4 receptor 
model A multiple sequence alignment was performed 
using the CLUSTALW algorithm [23]. A pi adrenergic 
receptor (A2A with a carvedilol agonist, PDB ID: 4AMJ) 
was used as the template because it had the highest se- 
quence identity percentage in the sequence alignment 
(41%). Protein models were generated from the alignment 
in a stepwise manner. The backbone coordinates for 
the aligned positions were extracted from the template 
and the regions of insertions/deletions in the alignment 
were found by searching either a loop library or a con- 
formational space search using constraint space program- 
ming. The templates were weighted by their sequence 
similarity to the target sequence and outlier atomic posi- 
tions were excluded. The scoring function used for asses- 
sing favorable interactions (hydrogen bonds, disulfide 
bridges) and unfavorable close contacts for determining 
side chain conformations was derived from a backbone- 
dependent rotamer library [24]. 



MD simulations 

All-atom MD simulations were carried out using CHARMM 
(v. c36a2) [25]. All simulations were carried out at 
323 K and 1 atm using periodic boundary conditions 
(PBC) with the NPT ensemble. The particle mesh 
Ewald (PME) algorithm was used for the long-range 
electrostatic interactions. Both the A2A receptor from 
PDB coordinates [14] and a model of the hERGl channel 
[1,2,26] were embedded into the DPPC membrane bilayer 
using the CHARMM-GUI membrane builder protocol 
[27]. Structures were minimized and equilibrated with 
gradually decreasing harmonic constraints (i.e., they 
were initialized with 10.0 and 5.0 kcal mol"^ A"^ for the 
backbone and side chains, respectively, and gradually 
decreased to 0.5 and 0.1 kcal moP^ A"^, respectively) 
over 2 ns (for equilibration) and then subjected to a 
50 ns production run. 

MM/GBSA computations 

The enthalpies of cisapride binding to both targets were 
computed as averages from an ensemble of structures 
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Figure 3 (left) Dynamics of cisapride in the hERGl cavity (Baclcbone RMSD after 50 ns ~5.0 A for a tetramer, 4.2 A for a monomer; the 
pore domain has an RMSD of ~2.9 A), (right) Summary of the MM/GB-SA decomposition analysis (per-residue and per-group contributions) 
from the last 40 ns of the MD simulation of the hERGl -cisapride complex, (see Figure 2 for a definition of the R groups at cisapride). 
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(5000 for each binding pose) sampled from evenly distrib- 
uted points over five independent 50 ns runs. The electro- 
static contributions to the desolvation components of the 
binding enthalpies were obtained by using the Generalized 
Born Solvation Energy Module and an Implicit Membrane 
(GBIM) module as implemented in CHARMM [28,29]. A 
dielectric constant of 2 was assigned to the protein, and 
the protein-solvent surface was defined using a set of opti- 
mized atomic radii from Nina et al [30]. Following the 
numerical recipe of Chandra et al [31], the membrane 
was represented as a 24 A slab with a dielectric constant 
of 2. To better describe the lipid dynamics, we extracted 
the positions of the lipids heavy atoms in the bilayer and 
used them to represent the neutral slab around the protein 
and cisapride. The protein-ligand interaction components 
of the binding enthalpies were computed following a 
standard protocol with infinite cut-offs [12,32]. 

Results and discussions 

MD simulations 

As stated in the Introduction, one of the major goals of 
this study was to identify the determinants of the 



cisapride blockade of hERG and its binding to the A2A 
receptor which serves as a suitable model for the 5HT-4 
receptor. Although our docking studies unambiguously 
identified binding pockets in both proteins, the orienta- 
tion of the bound ligand was less well defined. To ac- 
count for the limited resolution of our docking studies, 
we ran five independent 50 ns full-membrane all-atom 
molecular dynamics (MD) simulations using the top 5 
poses from docking analysis as a starting point. All sim- 
ulations were found to produce stable trajectories with 
membrane proteins displaying heavy atom RMSD values 
of ~ 3 to 4 A, which is comparable to root mean squared 
fluctuation (RMSF) estimates from previous studies of 
K-channels with available crystal structures [33]. The av- 
eraged locations of the ligand after MD simulations are 
shown in Figures 3 and 4. Cisapride bound in either 
membrane protein (A2A receptor or hERGl open-state 
channel model) displays significant conformational dy- 
namics in the binding site with an average RMSD value 
of -1.1 A relative to the average structure. The back- 
bone RMSD of the hERG tetramer channel after 50 ns 
was approximately 5.3 A including deviations from the 
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Figure 4 (left) The dynamics of cisapride in the A2A binding poclcet (Baclcbone RMSD after 50 ns ~3.1 A), (right) Summary of the MM/GB-SA 
decomposition analysis (per-residue and per-group contributions) from the last 40 ns of the MD simulation of the A2A-cisapride complex. 
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symmetric tetramer; the corresponding values for the 
monomers were from 3.5 to 4.2 A. The positional fluctu- 
ations in backbone atoms of residues forming the pore 
domain (PD) of the hERGl monomer plateaued at a 
RMSD of -2.4-2.9 A. These values are similar to those 
reported previously for MD simulations of Kv channels 
[32]. The backbone RMSDs of the A2A receptor plat- 
eaued at -3.1 A. The considerable dynamics of the re- 
ceptor and the bound ligand suggest that explicitly 
accounting for site and ligand flexibility in evaluation of 
the binding energies is necessary. To circumvent obvious 
limitations of the chosen docking strategy and to identify 
key interacting partners e.g. relevant amino-acid residues 
and key functional groups in the drug, we performed 
MM/GBSA computations to estimate the binding affinity 
of cisapride in the A2A receptor and the open-state hERGl 
model. The distribution of binding enthalpies from 5 
independent simulations is shown in the Figure 5. The 
average binding enthalpies for cisapride to the hERGl 
and A2A are -21.3 ± 2.8 kcal/mol and -24.3 ± 1.9 kcal/ 
mol, respectively. 

Analysis of the per-residue and per-functional group 
contribution to binding enthalpies allows us to deter- 
mine residues and drug fragments, which uniquely de- 
termines the binding mode of cisapride in both systems. 
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Figure 5 Normalized probability distribution function of the 
binding enthalpies for cisapride to WT-A2A and WT-hERGI proteins. 




(Figures 3 and 4) The key residues that determine the 
thermodynamics of cisapride stabilization in the hERGl 
cavity are S624 from the pore helix, Y652 and F656 from 
the distal S6 and, surprisingly, W568 from the S5 helix. 
The high-affinity intra-cavity blockade of hERGl by cisa- 
pride depends on the hydrogen bonds formed by the 
pore-helix and on stabilizing aromatic interactions with 
residues in the S5 and S6 helices. The drug dimensions 
allow for inter-subunit interactions (see Figure 2) with 
matching hydrophobic moieties in the hERGl cavity while 
the carbonyl backbone of cisapride interacts strongly 
with T623 and S624. The binding mode is similar to that 
for dofetiUde and other molecules known to restrict 
hERGl currents by binding to the proteins internal 
cavity [3,26]. Most of the residues were mapped in previ- 
ous experimental studies and are thought to be respon- 
sible for high-affinity blockades [34-38]. Importantly for 
both systems, docking studies were able to identify key 
amino acids and functional groups and thus results ob- 
tained from average protein structure and flexible ligand 
docking appear to be sufficiently predictive to allow for 
micromanagement of the drug. 

The decomposition of binding enthalpies for cisa- 
pride binding to the A2A adenosine receptor is iUus- 
trated in Figure 4. While there is a cluster of acidic 
residues involved in drug binding (El 3, El 69 and 
D170), the desolvation penalty for the negative charges 
located on the surface of the receptor counter balances 
favorable charge interactions and contributes unfavor- 
ably overall to drug stabilization in the pocket. Instead, 
the bound cisapride is stabilized to a large extent by a 
network of aromatic and amphipathic groups lining the 
binding pocket. The constellation of residues includes 
Y271, L267, 1252, T256 and M270 that interact with the 
linker and the rings (Rl, R3, and R6). (See Figure X for 
labeling of the regions of cisapride) The binding of cisa- 
pride is further stabflized by interactions between the 
residues S67, T68 and F168, and the moieties of the R2 
and R3 regions. It is worth mentioning that all of these 
residues are involved in the formation of binding 
pockets for adenosine and other agonists, according to 
previous studies and recently published crystal struc- 
tures [14]. F168 forms binding interactions with the 
central methoxy groups of the bound molecule, and 
other key residues (S67, T68 and Y271) interact with 
amphipathic functional groups on aromatic rings. The 
molecular fragments of the drug that form stable inter- 
actions with the residues in the A2A receptor are as fol- 
lows: a methoxy group on the central heterocyclic ring 
(R4), a terminal aromatic ring (Rl), and a methoxy 
group (R7) on the benzamide ring. Interestingly, it has 
been shown that modification of the benzamide ring is 
essential for modulation of the activity of cisapride de- 
rivatives [39]. 
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To better understand the effects of each fragment of 
cisapride on the A2A adenosine receptor binding and 
blocking affinities of hERGl, we also performed per-group 
decomposition of the binding enthalpies for drug binding 
to hERGl. The cisapride fragments were defined accord- 
ing to the topology definition of the CHARMM general 
force field [25,40] and are illustrated in Figure 2. The 
resulting decomposition is shown in Figures 3 and 4, 
allowing us to focus on targeted modifications of the frag- 
ments. It is apparent that to inhibit binding to hERGl, one 
may need to modify one or both of the terminal rings (Rl, 
R6) as well as the flexible linkers (R2, R3). Cisapride 
bound to the hERGl cavity interacts with a ring of hydro- 
phobic and aromatic residues (Y652, A653 and F656) and 
the Ser/Thr-rich apex of the pore helix (T623, S624) [2]. 
Experimental studies have also demonstrated that the 
matching interactions between cisapride and aromatic 
rings (F656 and Y652) are essential for the hERG blockade 
[41]. The decomposition analysis shows that one of the 
key interactions important for the hERG blockade but un- 
related to the high-affinity binding of cisapride to the A2A 
receptor is the van-der-Waals interactions between the 
long alkyl chain linker of the drug and a number of hydro- 
phobic residues in the hERG cavity. The docking/MD 



simulations strategy adopted in this study appears to 
have yielded a potential route to rehabilitate cisapride. 
In our next step, we attempt to remove fragments re- 
sponsible for drug stabilization in the inner-cavity of 
hERGl channel, while retaining its binding to the A2A 
receptor. Several parallel approaches were considered 
for de novo development of non-blocking cisapride col- 
lected in Table 1: 

(i) In the first approach, we simply targeted 
heteroatoms in the functional groups of cisapride 
molecule. They were exchanged with hydrogen (-H) 
to observe the effect of each of these fragments on 
the docking scores for both receptors. 

(ii) Next, the number of -CH2 groups of the linker was 
modified to produce truncated and extended 
versions of the linker. This affects the distance 
between the terminal aromatic ring and the 
six-membered heterocyclic ring piperidine, thus 
affecting torsional energy. The length of the alkyl 
chain also modulates the interactions between the 
aromatic moieties of the drug and the Y652/F656 
residues that are key components of the high- 
affinity blockade of hERGl. 



Table 1 Selected cisapride analogues from the de Novo drug design study 



Ligands 



2D structures 



A2A receptor 

docking score (kcal/mole) 



hERGl 

docking score (kcal/mole) 



Cisapride 



Cisapride_Frag_337 




-9.46 



-4.73 



Cisapride_Frag_1 82 



3.66 



-5.34 




The 10700 small organic molecules fragment database of Schrodinger was used to create new derivatives of cisapride (the place for replacement was previously 
found by a trajectory analysis of the MD simulations). Molecules were initially docked using virtual high-throughput screening (VHTS) with subsequent Glide/XP 
docking for the promising derivatives (-500 compounds). Finally, the inhibitory profiles of these derivatives at hERGl were tested using Induced Fit Docking. 
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(iii) In the last approach, the functional groups were 
removed based on the per-fragment decomposition, 
as shown in Figures 1 and 3. 

A library of 10,700 molecules from the small organic 
molecules fragment database of Maestro was used to 
create new analogues of cisapride guided by energy de- 
composition analysis as well as the per-residue inter- 
action energy analysis described above. The summary of 
the modification strategy is shown in Additional file 1: 
Tables SI to S2. All of the developed molecules were ini- 
tially docked using virtual high throughput screening 
(VHTS) followed by Glide extra precision (XP) docking 
for the derivatives showing low hERG blocking binding 



scores and high binding scores for the A2A adenosine re- 
ceptor (-500 compounds). 

General note on molecular docking and De novo drug 
design 

All of the derived cisapride analogues were docked to 
the A2A (PDB ID: 3QAK) active site using the approach 
described in the methods section. They were also docked 
to the pore domain of hERGl model, and the cisapride 
derivative results were compared to the original cisa- 
pride docking scores. (Figures 6 and 7, and Additional 
file 1: Table SI) Four different molecular docking pro- 
grams (AutoDock, Glide, FlexX, and GOLD) were used 
to study the binding interactions of cisapride and the 




Durdagi et a/, BMC Pharmacology and Toxicology 2014, 15:14 
http://www.bionnedcentral.conn/2050-65 1 1 /1 5/1 4 



Page 9 of 1 5 





w 

246 



V 

84 



A 

59 



s 

277 



i 
249 



F 

52 



M 

174 



c 

166 



A 
63 



Y 
271 



S 

67 



L 
167 



69 267 



T 

68 




Figure 7 Ligand Binding to A2A receptor. Top panel: The top-docking pose of cisapride in the human A2A adenosine receptor, (left) Binding 
interactions are zoomed and detailed (right). Bottom panel: A 2D-ligand interaction map for cisapride binding to the A2A receptor. 



cisapride derivatives in the A2A adenosine receptor and 
the hERGl central cavity. Since the FlexX and AutoDock 
docking programs both underestimated the interaction 
energies of cisapride and its derivatives at both targets 
(i.e., hERG and A2a)> GOLD and Glide/XP are generally 
preferred for the evaluation of binding affinities of 
known compounds. Our previously developed model 
of the trans-membrane domains, S1-S6, of open-state 
hERGl [1,26] was used to model the binding interac- 
tions of cisapride and its derivatives in the pore domain 
of the channel. Initially, cisapride was docked to A2A 
receptor as well as hERGl model, and these binding 
scores were used as threshold values for all four 
docking programs. Cisapride derivatives that displayed 



similar/higher binding scores (absolute values) in the 
A2A receptor and lower binding scores in the hERGl 
pore domain were considered for further analysis. For 
example, one of them (#11) labeled as Cisapride-Dll 
in the text had docking scores similar to original cisapride 
in the A2A binding pocket and considerably lower docking 
scores in the hERG PD, thus, it was considered for further 
rehabilitation. 

Effect of functional group substitution on drug binding to 
A2A receptor and hERGl channel 

All of the candidates with hydrogen substitutions by het- 
eroatoms (i.e., Cisapride-Dl and Cisapride-D8 in the 
Additional file 1: Table SI) had docking scores similar to 
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those of cisapride binding to the A2A receptor. However 
all of these derivatives retained high docking scores for 
binding to the hERG PD and therefore were discarded. 
Some of the new derivatives, such as Cisapride-D2 (an 
exchange of a -CI group with -H), led to a decrease in 
the docking score in the A2A receptor (i.e., the Glide/XP 
docking score decreased from -8.48 to -5.07 kcal/mole 
for Cisapride-D2), thus these derivatives were also dis- 
carded. (Additional file 1: Table SI) While the removal 
of some functional groups from cisapride did not sig- 
nificantly reduce the docking scores of compounds at 
A2A> their hERG PD bindings were still too high (i.e., 
Cisapride-D4, Cisapride-D5). 

Varying the number of -CH2 groups in the linker, R3 
between the terminal aromatic ring, Rl, and the central 
heterocycUc ring, R4, affected the docking scores of 
compounds in the hERG PD. (See Figure 2 for labeling) 
While decreasing the number of -CH2 groups from 
three (original cisapride) to one reduced the docking 
score of the derivative in the hERGl central cavity, (i.e., 
for Cisapride-Dll, GOLD docking score went from -9.42 
to -7.40 kcal/mole; the Glide/XP docking score went 
from -7.80 to -7.09 kcal/mole) increasing the number 
of -CH2 groups from 3 to 4 did not significantly affect 
the docking score of the derivative in the hERGl PD 
(i.e., for Cisapride-D9, the GOLD docking score went 
from -9.42 to -9.03 kcal/mole; the Glide/XP docking 
score went from -7.80 to -7.45 kcal/mole) (Additional 
file 1: Table SI). Therefor for further rehabilitation 
studies, only a decreased number of -CH2 groups in 
the cisapride backbone (i.e., Cisapride-Dll) was con- 
sidered. Derivatives of Cisapride-Dll were generated 
using the Combinatorial Library Enumeration and Screen- 
ing module of the Maestro molecular modeling package 
and the corresponding docking results were compared 
with that of original Cisapride-Dll (Additional file 1: 
Table S2). For this purpose, approximately 600 default 
fragments in Schrodinger s Enumerate module were used, 
and 8600 new cisapride derivatives were generated. After 
the derivative generation, the compounds were prepared 
(protonation states were determined at a pH of 7) with 
the LigPrep module of Maestro and energy minimization 
was performed for the structures using the PRCG energy 
minimization method. Structures were then docked onto 
the A2A adenosine receptor via the Glide High Through- 
put Virtual Screening method (Glide/HTVS). This was 
used to reduce the large number of derivatives to a man- 
ageable number for further analysis. The top 100 com- 
pounds according to the Glide/HTVS docking scores 
were used for Glide extra precision (XP) docking both in 
the A2A active site and the hERGl central cavity. The 
docking results of the selected derivatives of Cisapride- 
Dll are tabulated in the Additional file 1: Table S2. One 
of the major outcomes of this in silico study is that simple 



truncation of the flexible linker appears to be sufficient to 
remediate unwanted off-target interactions of cisapride. 
Obviously, it is not sufficient just to claim that such modi- 
fications will lead to a potent drug; it requires verification. 
One of the reasons for choosing cisapride is the availability 
of a large number of previously evaluated cisapride deriva- 
tives. Accordingly, to test similarities and differences in 
predicted cardiotoxicity (hERGl blockade) and efficacy of 
the cisapride derivatives we used the ZINC database and 
literature mining. 

Comparison study with a database of available 
drug databases 

94 cisapride derivatives (with >80% structural similarity 
with cisapride) taken from the ZINC databank were 
used to obtain the docking scores/poses for binding to 
the human adenosine A2A receptor and the hERGl 
pore domain using the Glide/XP docking program. The 
compounds that had docking scores (absolute values) 
greater than -7.00 kcal/mole in the A2A receptor bind- 
ing site (28 cisapride derivatives) were evaluated for 
their docking scores in the hERG PD (Additional file 1: 
Tables S3 and S4). Results showed that a few promising 
compounds (i.e., ZINC05998832, ZINC58529167, ZINC 
13834042, ZINC20621758, ZINC43023913) with docking 
scores less than (absolute values) -6.00 kcal/mole in the 
hERGl PD, and greater than (absolute values) -7.00 kcal/ 
mole in the binding pocket of the A2A receptor. (Table 2) 

Interestingly, ZINC20621758, also known commer- 
cially as mosapride, was detected via blinded dual-target 
screening to be a safer 5HT-4 agonist than cisapride. 
Potet et al [6] showed that mosapride does not block 
the hERG channel, which is in agreement with lower 
binding scores predicted in silico. It is also a well-known 
agonist of the 5HT-4 receptor. The IC50 of mosapride 
binding to the hERG channel is -16.5 (iM. Mosapride 
carries a significant decrease of cardiovascular risks re- 
lated to alterations in QT intervals according to preclin- 
ical studies [42,39,43]. Carlsson et al used a rabbit 
model of the acquired long QT syndrome, and while 
cisapride prolonged the QT interval, mosapride did not 
[39,44,45]. One of the key differences between mosa- 
pride and cisapride is a shorter linker between aromatic 
rings that, as predicted in the dual-target de novo drug 
modeling, will have an immediate effect on the drug- 
hERG PD interactions. Figure 8 shows the top Induced 
Fit docking pose of mosapride in the A2A receptor. Cisa- 
pride and mosapride have different binding modes in the 
A2A receptor due to a different number of CH2 groups 
in the linker. Additional file 1: Figure S2 shows a super- 
position of the top-docking poses of cisapride and 
mosapride. 

The other ZINC compounds in Table 2 also exhibit 
high docking scores in the A2A receptor; however there 
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Table 2 Cisapride derivatives from the ZINC database that show low binding scores in the hERG1 PD and high binding 
scores in the A2A adenosine receptor 



Cisapride derivatives 2D structure 



Glide XP docking Glide XP docking 

score (S)A2A (kcal/mol) score (S)hERG1 PD (kcal/mol) 



ZINC05998832 



ZINC58529167 






-9.45 



-7.57 



-4.62 



-4.51 



ZINC! 3834042 




-7.39 



-5.51 



ZINC20621758 




-7.22 



-5.67 



ZINC43023913 



-7.21 




-5.01 



are no experimental results available in the literature for 
these compounds. However, ZINC05998832 {aka fluoro- 
clebopride) is a close analog of clebopride. The reported 
IC50 of clebopride in the hERG PD is 0.62 (iM. Clebo- 
pride is commercially available in Spain and Italy as a 
gastroprokinetic drug. Although its cardiotoxicity has 
not been reported in clinical studies, in vitro studies 
have shown that clebopride prolongs the cardiac action 
potential duration at 90% (but not 50%) repolarization at 
10 (iM [44]. Tack et al [44] also found that no cardio- 
vascular safety concerns were reported for the newer se- 
lective 5HT-4 agonists prucalopride, velusetrag, and 
naronapride or for the non-selective 5HT-4 agonists 
with no hERG or 5HT-1 affinity, such as renzapride, cle- 
bopride, and mosapride. Importantly, the experimental 
data available for several ZINC compounds indicated nM 



affinity against 5HT-4 and 5HT-2 receptors, correspond- 
ing to high affinity binding for several of the compounds 
studied (ZINC20621758, Ki --113 nM; ZINC43023913, 79 
nM; ZINC13834042, 2250 nM and ZINC05998832 - 283 
nM, respectively), all with very low docking scores in the 
hERG PD. Although a few cisapride analogues (such as 
ZINC28087327, Ki 40.6 nM in 5HT-4; ZINC43024452, Ki 
3.16 nM in 5HT-4) were correctly predicted by in silico 
screening to have high experimental binding affinities in 
the 5HT-4 receptor, they also had high docking scores in 
the hERG PD (Additional file 1: Tables S2 and S3). 

On the transferability of A2A predictions to 5HT-4 receptors 

There is a well-known challenge in relating the top bind- 
ing poses produced by docking to the most stable ligand 
orientations in the binding pocket. Many organic 
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Figure 8 Mosapride bidning to hERGI pore domain. Top panel: The top-docking pose of mosapride that shows a low hERG PD blocking 
affinity in the human A2A adenosine receptor, (left) Binding interactions are zoomed and detailed (right). Bottom panel: 2D-ligand interaction 
map for mosapride binding to the A2A receptor. 



molecules are flexible with a number of rotating bonds. 
The previous section on mapping the relevant binding 
sites in the receptor suggests that some of these residues 
are also responsible for selective substrate and agonist 
binding to human 5HT-4 [15,46]. It is worth mentioning 
that the main target of cisapride in gastrointestinal tis- 
sues is not the A2A receptor but the 5HT-4 receptor. To 
investigate the similarities of the binding pockets of two 
GPCRs, namely, 5HT-4 and A2A receptors, we per- 
formed homology modeling with subsequent drug dock- 
ing. The homology modeling was performed using 
SWISS -MODEL for 5HT-4 (using the highest sequence 
identity scores of a multiple sequence alignment), pi ad- 
renergic receptors (with agonist bound) resulted in the 



highest sequence identity percentage from a sequence 
alignment (41%) and it was used as a template. Superpos- 
ition of the 5HT-4 model protein and the A2A crystal struc- 
ture reconfirmed that there is high structural similarity 
between these two GPCRs (RMSD of 2.4 A) (Additional 
file 1: Figure SI). Furthermore, most of the residues identi- 
fied by MD simulations as essential for drug stabilization 
are well-known determinants of high-affinity binding in 
the 5HT-4 system. We also used the ChEMBL database to 
search for compounds targeting the 5HT-4 receptor 
(CHEMBL1875). ChEMBL is a database for bioactive 
drug-like small molecules that contains abstracted experi- 
mental bioactivities from the scientific literature (i.e., bind- 
ing constants and pharmacology and ADME/Tox data) 
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[34,35]. In total, 419 associated bioactivity results for the 
5HT-4 receptor were found. The results were ranked based 
on Ki values, and compounds that had better Ki values 
than 1 nM in the 5HT-4 receptor were kept for consider- 
ation. Our aim was to test the docking score accuracies for 
compounds in the binding pockets of the A2A and 5HT-4 
receptors. The protonation states were determined at 
pH 7, and structures were optimized using the LigPrep 
module of Maestro. Subsequently, these compounds were 
docked to the A2A and 5HT-4 receptors using Glide/XP 
Induced Fit Docking. Considering the simplicity of the 
docking approach, our results show a high correlation 
between the experimental and predicted binding ener- 
gies (r^ > 0.6 for both compounds in the A2A and 5HT- 



4 receptors, docking results. Figure 9). Structures with 
bioactivity results in the 5HT-4 receptor were also 
docked in the hERG PD using Induced Fit Docking. 
Compounds with low predicted binding affinities in the 
hERG PD are listed in Additional file 1: Table S5 and 
can guide future studies. Our findings reconfirm the 
similarity of the binding pockets of the A2A and 5HT-4 
receptors as well as the transferability of our results. 

Conclusions 

A cross analysis of the interactions between cisapride 
and its analogues with the human A2A adenosine recep- 
tor and the hERGl central cavity led us to formulate a 
computational approach to the rehabilitation of drugs 



Co imparl sen of Experimental and Calculated Binding Energies of Compounds at A2A Receptor 





Comparison of EMperimental and Calculated Binding Energies of Compounds at 5HT4 




Figure 9 Comparison of experimental and calculated (IFD docking scores) binding energies of compounds. Top panel: The A2A receptor 
and Bottom panel: The 5HT-4 receptor. (Ki values are converted to kcal/mole for comparison, T = 300 K, R = 8.314 J/molK; If a compound has 
different Ki values for same target, the average of these values is used). 
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withdrawn from the market due to pro-arrhythmic activ- 
ity. Our SAR study of cisapride and its derivatives allows 
for the in silico evaluation of drug potency and the 
drugs ability to block hERG. Molecular docking results 
using both de novo designed compounds and available 
cisapride derivatives in the ZINC drug database showed 
that the shorter alkyl chains in the cisapride analogues is 
key element to retaining their binding to the A2A recep- 
tor and remediating the blockade of the hERGl channel. 
We believe that a simple dual-target 'rehabilitation 
strategy based on an optimization against two protein 
target structures may be applied to other drugs with- 
drawn from the market due to their side effects, and 
may lead to the reuse of these drugs. 

Additional file 



Additional file 1: Summary of Docking, De Novo Design and 
Modeling studies for cisapride analogues investigated. 
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